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ABSTRACT 

We examine the dependence of the mass-to-light (M/L) ratio of large-scale 
structure on cosmological parameters, in models that are constrained to match 
observations of the projected galaxy correlation function w p {r p ) and the galaxy 
luminosity function. For a sequence of cosmological models with a fixed, ob- 
servationally motivated power spectrum shape and increasing normalization erg, 
(SI ' we find parameters of the galaxy halo occupation distribution (HOD) that re- 



produce w p (r p ) measurements as a function of luminosity from the Sloan Digital 
Sky Survey (SDSS). From these HOD models we calculate the r-band conditional 
luminosity function $(L|M^), and from this the mean M/L ratio as a function 
of halo mass M^. We also use $(L|M/j) to populate halos of N-body simulations 
with galaxies and thereby compute Mj L in a range of large-scale environments, 
including cluster infall regions. For all cosmological models, the M/L ratio in 
high mass halos or high density regions is approximately independent of halo mass 
or smoothing scale. However, the "plateau" value of Mj L depends on <7g in addi- 
tion to the obvious proportionality with the matter density parameter Q m , and it 
represents the universal value (M/L) = VL m p CT - lX ,/ p\ um only for models in which the 
galaxy correlation function is approximately unbiased, i.e., with <7g ~ a% g . Our 
results for cluster mass halos follow the trend (M/L) cl = 577 (f2 m /0.3) (crg/0.9) 17 
/iM /L . Combined with Carlberg et al.'s (1996) mean M/L ratio for CNOC 
galaxy clusters, this relation implies (a 8 /0.9) (fi m /0.3) 06 = 0.75 ± 0.06. M/L 
estimates for SDSS clusters and the virial regions of clusters in the CAIRNS 
survey imply a similar value of <7gf2^; 6 , while the CAIRNS estimates for cluster 
infall regions imply a lower value. These results are inconsistent with param- 
eter values Q m ~ 0.3, erg ~ 0.9 favored by recent joint analyses of cosmic mi- 
crowave background measurements and other large scale structure data, though 
they agree with values inferred from van den Bosch et al.'s (2003) analysis of the 
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2dF Galaxy Redshift Survey. We discuss possible resolutions of this discrepancy, 
none of which seems entirely satisfactory. In appendices we present an improved 
formula for halo bias factors calibrated on our 360 3 N-body simulations and an 
improved analytic technique for calculating the galaxy correlation function from 
a given cosmo logical model and HOD. 

Subject headings: cosmology: observations — cosmology: theory — galaxies: 
clustering — galaxies: clusters — large-scale structure of universe 

1. Introduction 

Determining the matter density of the universe is one of the key goals of observational 
cosmology. By definition, the mean matter density is the product of the mean luminos- 
ity density and the mean mass-to-light ratio of the universe, (M/L), making the density 
parameter 

Qm = (M/L) X pi um /pcrit- (1) 

One of the classic methods of inferring VL m is to estimate the Mj L ratios of galaxy groups or 
clusters, then multiply by the mean luminosity density derived from the galaxy luminosity 
function (e.g., Gott et al. 1974; Peebles 1986; Bahcall et al. 1995; Carlberg et al. 1996). 
However, dynamical methods of estimating M/L necessarily focus on dense regions of the 
galaxy distribution, so this route to Q m relies on the assumption that the galaxy population in 
these region is representative of the universe as a whole. Observed M/L ratios rise steadily 
from the scale of binary galaxies to groups to rich clusters, but there is an approximate 
plateau in M/L values for the richest bound systems, and the mass-to-light ratios inferred 
from the dynamics of superclusters and the extended infall regions around clusters extend 
this plateau to larger scales (see, e.g., Bahcall et al. 1995; Bahcall et al. 2000, hereafter BOO; 
Rines et al. 2004). The existence of this plateau is sometimes taken as evidence that the 
measured values of M/L do indeed represent the universal value. 

This paper has two goals. The first is to assess the above line of reasoning: does the 
existence of a plateau in M/ L at large scales imply that the universal value has been reached? 
The second is to assess the implications of observed M/L ratios for the values of Vt m and 
the matter fluctuation amplitude cr$, which we will show plays a crucial role. (Here a% is 
the rms matter fluctuation in spheres of radius 8 h" 1 Mpc, calculated from the linear matter 
power spectrum, where h = H /100 kms -1 Mpc -1 .) The key to these assessments is a model 
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for the relation between galaxies and dark matter (a.k.a. galaxy bias) that extends to the 
non-linear regime. We will derive this relation empirically, by fitting models of the halo 
occupation distribution (HOD; see Berlind & Weinberg 2002 and references therein) to the 
projected correlation functions of galaxies in the Sloan Digital Sky Survey (SDSS; York et 
al. 2000), along the lines of Zehavi et al. (2004a,b). HOD models describe bias at the level of 
virialized halos by the probability P(N\M) that a halo of virial mass M contains iV galaxies 
of a specified type, together with prescriptions for spatial and velocity biases within halos. 

BOO also addressed the complications of mass-to-light ratios, drawing on the results of 
a hydrodynamic simulation of the galaxy population. Their simulation predicts a plateau 
in the 5-band M/L that is close to the universal value, except for a modest "anti-bias" 
arising from the older stellar populations of galaxies in dense environments. Applying this 
anti-bias correction to the observed (M/L) B values, they infer fl m = 0.16 ± 0.03. However, 
BOO considered only a single cosmological model, with a 8 = 0.8. We will show that, when 
the bias relation is inferred empirically by matching the observed correlation function, a 
plateau in M/L at high halo masses or overdensities is a generic result. However, the 
plateau corresponds to the true universal M/L only for models with er 8 m <7g g , where 0Sg is 
the rms galaxy count fluctuation. For lower <jg, galaxies must be positively biased to match 
the observed clustering, and the M/L plateau is below the universal value. For higher <j 8 , 
the galaxy population is anti-biased in the dense regions, and the M/L plateau is higher 
than the universal value. 

Zehavi et al. (2004b) infer HOD parameters as a function of galaxy r-band luminosity 
by fitting measurements of the projected galaxy correlation function w p (r p ), assuming a 
cosmological model with <r 8 = 0.9. We carry out a similar procedure, but we consider a range 
of cr 8 values from 0.6 to 0.95. We keep the shape of the linear matter power spectrum P\ in (k) 
fixed, regardless of ag or Q m , because it is well constrained empirically by the combination 
of microwave background anisotropies and large-scale galaxy clustering measurements (e.g., 
Percival et al. 2002; Spergel et al. 2003; Tegmark et al. 2004), and because the physical effects 
of changes in fl m are much more transparent if they are kept separate from Pn n (k) changes 1 . 
For a specified fl m and <Tg, we can calculate M/L as a function of halo mass analytically, 
given the HOD parameters derived by fitting w p {r p ) (see §2). To compute M/L over larger 
scales, we use the derived HODs to populate the halos of N-body simulations. 

As shown by Zheng et al. (2002), for a fixed linear power spectrum, a change to Vt rn 



1 In the inflationary cold dark matter (CDM) framework, a change in TO with all other parameters held 
constant alters the shape of Pu n (k). However, variations other parameters, such as the Hubble constant, the 
inflationary spectral index, spectral running, the baryon density, and the neutrino mass can compensate for 
those changes, at least to some degree. 
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simply shifts the mass and velocity scales of virialized halos, with negligible effect on the 
shape of the halo mass function or the clustering of halos at a given (scaled) mass. For two 
cosmo logical models that differ only in the value of Q m , we can obtain identical real-space 
galaxy clustering by simply shifting the HOD mass scale in proportion to Q m . With bias 
constrained by the observed clustering, therefore, our predicted M/L ratios for a given erg 
are simply proportional to Q m . We will often indicate this scaling by quoting predicted 
M/L ratios in units of uj . 3 h M Q /L , where u , 3 = f2 m /0.3. The scaling with h arises in 
the observations because inferred dynamical masses scale as h' 1 and luminosities as h~ 2 . 
The same h scaling arises in the predictions because the characteristic mass in the halo mass 
function is proportional to h~ l at fixed Q m and cr 8 and we automatically impose the observed 
luminosity scale, which is proportional to h~ 2 . Uncertainties in the value of h therefore have 
no impact on our conclusions. 

Most recent applications of the mass-to-light method yield Vl m ks 0.15 — 0.20 (e.g., 
Bahcall et al. 1995; Carlberg et al. 1996; BOO; Lin et al. 2003; Rines et al. 2004). Alternative 
methods of inferring Q m , from cosmic microwave background (CMB) anisotropies, the shape 
of the galaxy power spectrum, Type la supernovae, and the baryon-to-total mass ratio in 
clusters, have begun to converge on a significantly higher value of Q m ~ 0.3 (e.g. Turner 
2002). Recently Tegmark et al. (2004) combined WMAP CMB data with the SDSS galaxy 
power spectrum to infer Vl m = 0.30 ± 0.04. Seljak et al. (2004) combined these data with 
the SDSS Lyman-a forest spectrum to obtain Vl m = 0.28 ± 0.02. The tension between M/L 
estimates of Vt m and the higher values from other methods warrants a careful investigation 
of the assumptions that underly the M/L approach. We will show that the conflict can be 
alleviated if as is low, ~ 0.6 — 0.7. However, the discrepancy persists for a matter fluctuation 
amplitude of a 8 « 0.9, which is favored by the Tegmark et al. (2004) and Seljak et al. (2004) 
analyses. 

Our approach is similar in spirit to the conditional luminosity function (CLF) analyses 
of clustering in the Two-Degree Field Galaxy Redshift Survey (2dFGRS; Colless et al. 2001) 
by Yang et al. (2003; 2004), van den Bosch et al. (2003), and Mo et al. (2004). In practice 
there are important differences in the two techniques. In the CLF, the luminosity function 
within an isolated halo is parameterized by a Schechter function, while in our approach this 
function is determined non-parametrically. We fit the full projected correlation function, 
while the above CLF papers fit the luminosity dependence of the galaxy correlation length, 
together with the galaxy luminosity function. Despite the many differences in analysis and 
parameterization, and the use of independent galaxy clustering measurements derived from 
an r-band selected sample instead of a fej-selected sample, we reach a bottom line similar 
to that of van den Bosch et al. (2003): obtaining cluster mass-to-light ratios close to the 
observational estimates of ~ 350 hM Q /L Q (e.g. Carlberg et al. 1996; BOO) requires either 
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Vt m or <t 8 to be significantly below the currently favored values of 0.3 and 0.9 (see §4). 



The key ingredient in the HOD prescription of galaxy bias is P(N\M), the probability 
that a halo of mass M contains N galaxies of a specified class. The galaxy classes we will 
consider are defined by thresholds in r-band luminosity. Our model for P(N\M) is motivated 
by theoretical studies of the HOD using semi-analytic galaxy formation methods, high- 
resolution N-body simulations, and full hydrodynamic cosmological simulations (specifically 
Kravtsov et al. 2004 and Zheng et al. 2004; for earlier work see Jing et al. 1998; Kauffmann 
et al. 1997, 1999; Benson et al. 2000; Seljak 2000; White et al. 2001; Yoshikawa et al. 
2001; Berlind et al. 2003). This model distinguishes halo central galaxies from satellites in 
halos containing multiple galaxies. In our standard parameterization, the number of central 
galaxies above the luminosity threshold changes sharply from zero to one at a minimum halo 
mass M min . We also consider a model in which the probability of hosting a central galaxy is 
exp(— M min /M), to test the sensitivity of our results to the assumed sharpness of the central 
galaxy threshold. 

We adopt a functional form 



for the mean number of satellites in a halo of mass M > M min , and (N sat )M = for 
M < M min . Here M cut is a cutoff mass scale for the satellite galaxy power law, which allows 
a soft transition to halo masses that host no satellites. Zehavi et al. (2004b) showed that a 
model of this form allows good fits to the observed w p (r p ) for the concordance cosmology, 
though they mainly focused on an alternative parameterization with a sharp cutoff in (N sat ) M 
and a free power-law index (N sat ) M oc M a . The fixed-a with varying M cut parameterization 
is better suited to our purposes here because we will difference occupation functions for 
adjacent luminosity thresholds to obtain P(N\M) for luminosity bins, and small statistical 
errors in a can drastically affect the number of of satellites in individual bins at high halo 
mass. We assume a Poisson distribution of satellite number relative to the mean {N saX )M, 
as implied by the theoretical studies. 

This parameterization of P(N\M), while restricted, has the virtue of capturing theo- 
retical predictions quite well while introducing only three free parameters. An essentially 
perfect fit to theoretical results requires two additional parameters providing freedom in a 



2. HOD Modeling of w p {r p ) 




(2) 
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and the central galaxy cutoff shape (Zheng et al. 2004). However, we expect the three pa- 
rameter model to be adequate to our purposes here, predicting the mass-to-light ratios in 
and around high mass halos. 

We will investigate the effects of changing the assumed value of a, but we note that 
the choice of a = 1 is well motivated by the results of semi-analytic models, hydrodynamic 
simulations, and collisionless N-body simulations (Zheng et al. 2004; Kravtsov et al. 2004). 
Observational estimates of a are subject to the effects of interlopers, completeness limits, and 
uncertainty in the halo mass estimates, and the agreement among observational estimates is 
not as good as for the theoretical studies. Lin et al. (2004), using a sample of rich clusters 
from the 2MASS survey, find N sat oc M 87 , in agreement with the much smaller sample of 
Rines et al. (2004), who find a = 0.84. Kochanek et al. (2003), who analyze a different 
2MASS cluster sample, find a = 1.1. Analysis of the 2dF Percolated- Inferred Galaxy Group 
catalog (Eke et al. 2004) by Collister & Lahav (2004), which constitutes the largest group 
sample analyzed to date, yields a = 0.99, a value consistent with the preliminary results of 
the group multiplicity function for SDSS galaxies (A. Berlind, private communication). 

We fit the observed w p (r p ) for each luminosity threshold sample of SDSS galaxies mea- 
sured in Zehavi et al. (2004b), as listed in their Table 2. For the M r = —20 threshold 2 , we 
use the sample limited to redshift z < 0.06. We fit each sample for five values of a%: 0.95, 
0.9, 0.8, 0.7 and 0.6. Our results for erg — 0.9 are similar to those of Zehavi et al. (2004b), 
except for the slight differences in parameterization and calculational method (see Appendix 
B). For all calculations, we adopt a linear theory matter power spectrum with inflationary 
spectral index n s — 1 and shape parameter V = 0.2 in the parameterization of Efstathiou 
et al. (1992). Satellite galaxies are assumed to follow the "universal" halo density profile of 
Navarro, Frenk, & White (1997) with no internal spatial bias. The halo concentrations are 
calculated in the same manner as Kuhlen et al. (2004) assuming Vl m = 0.3, with the halo 
edge defined as the radius at which the mean interior density is 200 times the mean value. 
The calculation of w p (r p ) for a specified cosmology and HOD is described in Appendix B. 
Our method is similar to that of Zheng (2004), which was used by Zehavi et al. (2004b), but 
we implement an improved treatment of halo exclusion and adopt more accurate halo bias 
factors inferred from our N-body simulations (Appendix A). 

The data for each luminosity threshold are fit by minimizing x 2 using the full covariance 
matrix for each sample. The jackknife estimates of these matrices are discussed in detail 
by Zehavi et al. (2004b). The number of free parameters in the fit is reduced from three to 
two by matching the space density of galaxies for each sample. In practice, this constraint 



2 Throughout this paper, we quote absolute magnitudes for h = 1.0. 
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is used to fix the value of M min for a given M cut and Mi. We ignore the uncertainty in the 
space density itself. We require that the number density of satellite galaxies between two 
adjacent magnitude limits be smaller for the brighter sample at all halo masses, a consistency 
condition that only affected the result in two of the 45 fits performed. The results of the fits 
are listed in Table I. 

Figure la shows two examples of these fits. The observed w p (r p ) for M r < —20 and 
M r < —21.5 are shown by the symbols with error bars, plotted on top of the HOD fits 
for each value of cr 8 . For M r < —21.5, the different fits are difficult to distinguish. For 
M r < —20, the larger error bars allow the best-fit HOD models to differentiate some at 
large scales, even though the x 2 values per degree of freedom are all less than one. Note 
that strong covariances between data points make a simple visual estimate of x 2 unreliable. 
In Figure lb, the data are plotted against the nonlinear matter correlation functions for 
each value of erg. These correlation functions are calculated by Fourier transformation of 
the Smith et al. (2003) non-linear P(k), which takes the linear power spectrum as input, 
and assume Vt m = 0.3. The differences in the matter distributions are large, but are easily 
overcome even with our restrictive three-parameter HOD. 

The values of x 2 per degree of freedom for the HOD fits fluctuate from sample to sample, 
with a median value below 1.0 for <r 8 > 0.8 and above 1.0 for <r 8 < 0.7. We do not regard 
the somewhat higher x 2 values as significant arguments against low og, however, because 
our three-parameter description of the HOD is quite restrictive. In addition, adopting a 
CMBFAST (Seljak & Zaldarriaga 1996) power spectrum with the Tegmark et al. (2004) 
cosmological parameters in place of the T = 0.2, Efstathiou et al. (1992) power spectrum has 
a noticeable effect on x 2 values, though it makes almost no difference to the best-fit HOD 
parameters themselves. Adding freedom to the HOD parameterization or changing the input 
power spectrum could thus have significant effect on the quality of the w p (r p ) fits, but we 
expect only minor influence on our mass-to-light ratio predictions because the overall shape 
of (N)m is well constrained by the number density and w p (r p ) measurements. 

From the HOD fits listed in Table 1, we can calculate a discrete estimate of the condi- 
tional luminosity function, ®(L\M h ), the luminosity function of galaxies in halos of mass M h . 
Since the HODs describe samples brighter than a specified absolute magnitude, differencing 
the number of galaxies in a halo of mass Mh within the magnitude bin of width AM r = 0.5 
mag yields §(L\Mh) within that bin: 

HL\M h )AL = (JV)2J> - (N)^~ AM ^. (3) 
We use luminosity rather than magnitude on the left hand side of equation (3) for clarity 
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of notation. In this form, equation (3) is normalized such that a summation of $(L|M^) 
over all magnitude bins returns {N)m for the M r < —18 sample. The parameters in Table 1 
assume our standard HOD model with a sharp threshold for the central galaxy occupation. 
When we use an exponential cutoff for the central occupation, we keep the same values of 
Mi and M cut but adjust M min slightly so that the mean galaxy space density remains fixed. 

Figure 2 plots the conditional luminosity functions for three halo masses for the expo- 
nential central cutoff. Panel (a) shows <J>(L|M/j) for = 3 x 10 12 h^ 1 M Q , normalized so 
that the area under each curve is the total number of galaxies expected at this halo mass. 
The dashed line that falls rapidly at M r < —19.5 is the satellite galaxy contribution to 
§(L\Mh) for the model with as = 0.95. At this low halo mass, the number of satellite galax- 
ies brighter than M r = —18 is smaller than the number of central galaxies, but satellites 
dominate at the faintest magnitudes. The variation of $(L|M fe ) with <r 8 is minimal because 
at this point on the halo mass function the space density of halos is only weakly sensitive to 
<7g. The difference in the total number of galaxies brighter than M r = — 18 between the low 
and high values of a 8 is 1.84 — 1.56 = 0.28. 

Panel (b) shows §(L\M h ) for M h = 3 x 10 13 h~ l M . Halos of this mass are intermediate 
between those that host galaxy groups and clusters, with (N) M ~ 8 for M r < — 18. Satellite 
galaxies dominate the conditional luminosity functions out to M r = —21, beyond which the 
brighter central galaxies dominate. The different values of <jg now produce small differences 
in the overall normalization of $(L|M^). As a 8 decreases, the space density and clustering 
amplitude of high mass halos decreases, and the mean occupation of these halos must grow 
to keep the galaxy number density and clustering amplitude constant. In panel (c), which 
plots §(L\M h ) for cluster-sized halos of M h = 3 x 10 14 h~ l M , the differences in the models 
are clear. The low-(T 8 model has nearly twice as many galaxies per halo as a 8 = 0.95. For all 
the models, §(L\M h ) is dominated by satellite galaxies at all but the brightest magnitude 
bin. 

Our eventual qualitative conclusion is already evident from Figure 2. With the galaxy 
space density and clustering amplitude fixed to match observations, a low-a 8 model must 
have a larger fraction of its galaxies in massive halos, and it therefore predicts lower mass- 
to-light ratios in these halos. 

To quantify this point, we calculate mass-to-light ratios as a function of halo mass by 
summing over luminosity, 
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where % denotes the magnitude threshold, which runs from —18 to —22. All galaxies in each 
magnitude bin are assumed to be at the midpoint of the bin, and galaxies in the M r . < —22 
sample are assumed to have M r = —22.25, an assumption that has negligible influence on 
the results. Equation (4) is equivalent to M/L = M^j /^ max L $(L|M^) dL given our discrete 
estimate of $(L|M^) in equation (3). 

Before proceeding to our main results, we want to check the sensitivity of our calculations 
to the form of the central galaxy cutoff and the value of a. In Figure 3a, points with error 
bars are the SDSS data for M r < —20 from Zehavi et al. (2004b), and the solid line shows 
the fit with a hard cutoff and the parameters in Table 1 for erg = 0.95. Filled circles show 
Wp(r p ) for the same values of Mi and M cut , but an exponential cutoff, with M min adjusted 
to retain the galaxy density. The differences in the calculated correlation functions are 
barely discernible, mostly confined to the two-halo term at small scales, where the one-halo 
term dominates w p (r p ) anyway. Figure 3b shows an example of $(L|M/j) for the soft and 
hard central cutoffs for cr 8 = 0.95 at Mh = 3 x 10 13 h~ 1 M Q . The low luminosity regimes 
are identical, but at high luminosities where central galaxies dominate, the central cutoff 
produces a well-defined bump, which is smoothed out in the exponential cutoff model. Panel 
(c) shows the mass-to-light ratio as a function of halo mass for these two prescriptions. The 
sharp cutoff model produces artificial jumps in the M/L function, but the overall trend 
and, more importantly, the behavior at high halo masses are the same for hard and soft 
cutoffs. We conclude that the details of the central galaxy cutoff do not significantly affect 
the mass-to-light ratio predictions. 

In Figure 4 we quantify the dependence of our results on the value of a. First, we refit 
w p (r p ) and calculate M/L for cr 8 = 0.9 with a = 0.9, 0.95, 1.05, and 1.1. The results are 
shown, along with the fiducial results for a = 1, as the five curves in Figure 4a. The M/L 
ratios are nearly identical at group masses and below. At M 2 x 10 14 h~ l M the curves 
begin to separate, with lower a resulting in higher M/L, and vice versa, due to the different 
scalings of galaxy number with halo mass. 

Our second method for testing the sensitivity to a uses a much more flexible HOD 
parameterization in which the value of satellite mean occupation is specified at logM = 12, 
13, 14, and 15, connected by a cubic spline and smoothly truncated at low masses (see 
Figure 19 of Zehavi et al. 2004b for examples of similar fits). With this extra freedom, it 
is difficult to guarantee that (N S3 x)m is a monotonically increasing function of luminosity 
threshold, so we have therefore fit the Zehavi et al. (2004b) measurements for samples in 
absolute magnitude bins, rather than magnitude thresholds. In each bin, we assume that 
(N cen )M = 1 in a range M CiIllin — M Cjmax and (N ccu )m = elsewhere, equivalent to our sharp 
cutoff assumption for luminosity threshold samples. Using a Monte Carlo Markov Chain 
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(e.g., Lewis & Bridle 2002), we identify between 17 and 42 models for each luminosity bin 
that have Ax 2 < 1 relative to the best fit model. Using all 4.7 x 10 6 possible combinations 
of these models, we evaluate the M/L ratios at the cluster mass scale. The distribution 
of M/L values with respect to the mean at a given mass has an approximately Gaussian 
core with a half- width of 11% at 10% of the maximum. At lower probabilities, there is a 
tail toward high M/L values that arises from HOD fits that are statistically acceptable but 
physically implausible, with anomalously high halo masses for central galaxies. There is no 
corresponding tail to low M/L values, since if one puts too many galaxies in clusters the 
correlation function is inevitably too high. The shaded region in Figure 4a shows the full- 
width at 10% maximum as a function of halo mass. This test also accounts for statistical 
uncertainty in the w p (r p ) measurements themselves. 

Figure 4b quantifies the change in the average cluster M/L as a function of a. The y- 
axis shows the percentage difference in (M/L) d , averaged over the halo masses of the CNOC 
cluster sample (see §3.2 below), with respect to a — 1.0. For Aa = 0.1, A(M/L) C \ m 8%, 
and for Aa = 0.05, A(M/L) C \ ps 4%. Although the differences in the M/L curves of panel 
(a) appear large at the highest halo mass, the impact of these differences in moderated in 
panel (b) because clusters above 10 15 h^ 1 M are relatively rare. We conclude that for any 
given cosmology, the existing w p {r p ) measurements constrain the average cluster mass-to- 
light ratio, the quantity most relevant to our conclusions, to within about 10%. However, 
the detailed shape of the M/L curve at M h > 3 x 10 13 h~ l M depends on the assumed form 
of the HOD. 



3. Mass-to-Light Ratios of Large-Scale Structure 

3.1. M/L vs. overdensity 

The numerical results in this and subsequent sections are calculated using Tinker, Wein- 
berg, & Zheng's (2004) ACDM simulations (inflationary cold dark matter with a cosmological 
constant). These simulations, consisting of a set of 5 realizations of 360 3 particles evolved in 
a periodic volume 253 h^ 1 Mpc on a side, have the appropriate mass resolution for modeling 
the luminosity-dependent £ s (r) of galaxies with M r < —20. Fainter samples have minimum 
halo masses below the smallest resolved halos in our simulations, which are identified us- 
ing a friends-of-friends algorithm with a minimum of 30 particles. The dark matter power 
spectrum used to create the initial conditions is identical to that assumed in the analytic 
HOD calculations of the previous section. These simulations are run with gadget (Springel, 
Yoshida, & White 2000), with a force softening of 70 h^ 1 kpc, and the standard gadget 
timestep criterion rj = 0.2. The simulations have Q m = 0.1, A = 0.9, and a 8 = 0.95 at the 
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Table 1. HOD Parameters for the SDSS Galaxy Samples 



M r 


M min 


M cut 


Mi 


Xdof 


M min 


M cut 


Mi 


Xdof 


M min 


M cut 


Mi 


Xdof 






a 8 = 


0.95 






08 = 


=0.9 






a 8 = 


=0.8 




-18.0 


11.2 


12.0 


12.7 


1.63 


11.3 


11.7 


12.7 


1.40 


11.3 


12.1 


12.6 


l.n 


-18.5 


11.4 


12.3 


12.8 


0.88 


11.4 


12.3 


12.7 


0.72 


11.4 


12.7 


12.6 


0.49 


-19.0 


11.5 


12.7 


12.8 


0.94 


11.5 


12.7 


12.8 


0.90 


11.5 


12.9 


12.7 


1.04 


-19.5 


11.7 


13.1 


13.0 


0.22 


11.7 


13.1 


12.9 


0.22 


11.7 


13.1 


12.9 


0.51 


-20.0 


12.0 


13.0 


13.3 


0.51 


12.0 


13.0 


13.2 


0.51 


12.0 


13.2 


13.1 


0.58 


-20.5 


12.3 


13.1 


13.5 


1.80 


12.3 


13.2 


13.4 


2.21 


12.3 


13.3 


13.4 


3.09 


-21.0 


12.7 


13.6 


13.9 


1.28 


12.7 


13.6 


13.8 


1.36 


12.7 


13.7 


13.7 


1.58 


-21.5 


13.3 


14.6 


14.1 


0.94 


13.3 


14.6 


14.0 


0.90 


13.2 


14.6 


13.8 


0.86 


-22.0 


13.9 


14.8 


14.6 


1.42 


13.9 


14.7 


14.5 


1.42 


13.8 


14.9 


14.2 


1.27 






cr 8 = 


-fl 7 
-U. ( 






cr 8 = 


u.o 












-18.0 


11.3 


12.1 


12.5 


0.96 


11.3 


12.1 


12.4 


0.94 










-18.5 


11.4 


12.6 


12.6 


0.43 


11.4 


12.5 


12.5 


0.61 










-19.0 


11.5 


12.8 


12.7 


1.43 


11.5 


12.7 


12.6 


2.11 










-19.5 


11.7 


13.0 


12.8 


1.17 


11.7 


12.9 


12.7 


1.98 










-20.0 


12.0 


13.1 


13.0 


0.70 


12.0 


13.2 


12.9 


0.91 










-20.5 


12.3 


13.3 


13.3 


3.97 


12.2 


13.3 


13.1 


5.02 










-21.0 


12.7 


13.6 


13.6 


1.86 


12.6 


13.7 


13.4 


2.31 










-21.5 


13.2 


14.5 


13.6 


0.88 


13.1 


14.5 


13.4 


1.00 










-22.0 


13.7 


14.9 


13.8 


0.97 


13.6 


14.6 


13.8 


0.74 











Note. — Values of M min , M cut , and Mi are listed as \og(M/h~ l M Q ) 
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last output, and we use earlier outputs to represent z = results for lower ag, rescaling halo 
masses as necessary to represent different fl m (see Zheng et al. 2002 and Tinker et al. 2005 
for further discussion). For the purposes of this paper, this rescaling approach should yield 
nearly identical results to running separate simulations for each (fl m , a 8 ) combination, since 
at fixed cr 8 and Pn n (k) the value of Q m slightly alters the density profiles of halos but has 
minimal impact of the halo masses and clustering. 

Since $(L|M^) is known from the fits to w p (r p ), it is a simple process to create galaxy 
distributions from the dark matter halo populations in the simulations. These galaxy distri- 
butions will match both the SDSS clustering and luminosity function for all values of a$. For 
each halo above M min (for the M r < —20 sample), the central galaxy is placed at the center 
of the halo, which we define as the location of the most-bound dark matter particle, i.e., 
the particle with the most negative potential energy calculated by a Newtonian summation 
over all the particles in the halo. The number of satellite galaxies is chosen from a Poisson 
distribution with a mean given by equation (2), and the satellites are spatially distributed 
by random sampling of the dark matter particles in the halo (other than the most-bound 
particle). The luminosities of the galaxies are selected at random from $(L|M^), truncated 
at M r = —20. Since central and satellite galaxies have distinct contributions to $(L|M/j), 
in practice they are treated separately when populating the N-body halos; central galaxy 
luminosities are chosen only from their contribution to the total Q(L\M h ), while satellite 
luminosities are chosen from the satellite galaxy portion of <&(L\M h ). 

The halo mass function in the simulations differs from the analytic form of Jenkins et 
al. (2001) at some masses. The discrepancy is within the range shown in Jenkins et al.'s 
Figure 8, but at the high end, with maximum difference of 10-15% in number of halos at 
fixed mass. Using our analytically derived HOD parameters, which assume the Jenkins et 
al. (2001) form, therefore leads to ~ 7.5% too many galaxies (relative to the observed space 
density we are trying to match) when we populate the simulations. We correct for this 
discrepancy when populating the N-body simulations by changing M m i n for each sample to 
the value that yields the correct number of galaxies. This change is moderate, < 10% in 
^min, but it is required in order for the luminosity function in the simulations to accurately 
represent that of the Zehavi et al. (2004b) samples. 

A simple application of these galaxy distributions is to investigate mass-to-light ratios 
as a function of overdensity in randomly placed spheres, as done by BOO using their hydro- 
dynamic simulation. Figure 5 plots the mass-to-light ratio, relative to the mean M/L of the 
box, as a function of dark matter density at three top-hat smoothing scales. (In this way of 
plotting our data, values greater than one would be considered "anti-biased" by the defini- 
tion of BOO, since luminosity enters in the denominator). For each 5, the average light and 
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mass are calculated separately and then used to calculate M/L. This prevents divergence 
in spheres with no galaxies, but in practice it does not noticeably change the curves. Panel 
(a) shows the results for a top-hat smoothing scale of 1.5 h^ 1 Mpc, the scale used by Abell 
(1958) to define galaxy clusters, which is close to the virial radius of ~ 10 15 h" 1 M halos. 
The curves shoot up at 5 < 0, since at these overdensities there are no halos above M m \ n 
and therefore no galaxies, so p\ nm approaches zero faster than p m . The exact value of 5 at 
which this upturn occurs depends on the adopted luminosity threshold, since fainter galaxies 
occupy lower mass halos, and on the smoothing radius, since the probability of finding no 
halos above M min drops with increasing radius. The sharpness of the upturn depends on 
the sharpness of the {N) M threshold at M min , but the qualitative behavior at 5 < is a 
generic prediction of any model with a minimum halo mass. In the range < 5 < 10, all the 
curves dip below the mean due to low-mass halos that host a central galaxy. Overdensities 
of 5 > 100 begin to represent virialized, cluster-like structures. At this 5, the different a$ 
models spread out, with high-og having M/L ratios above or equal to the mean and \ow-a 8 
models having M/L as low as half the global mean. Panels (b) and (c) show the same 
results but for smoothing scales of 5 and 10 h^ 1 Mpc. For both of these smoothing scales, 
the different a 8 models diverge at the highest overdensities, with low a 8 corresponding to 
M/L ratios below the mean of simulation. At the highest overdensities, most of the M/L 
curves flatten to a roughly constant value, though for low a$ and large smoothing scale the 
M/L ratio is a declining function of 5 in this regime. More importantly, the values of M/L 
at high 5 are nearly independent of the smoothing scale for a given model, even when they 
are far from the universal value. 

The implications of Figure 5 are clear. The existence of a plateau in M/L, either as 
a function of 5 at fixed smoothing length or as a function of smoothing scale at fixed S, 
cannot be taken as evidence that one has measured the universal value of M/L. Simply 
measuring M/L in a large volume does not guarantee convergence to the universal {M/L) 
if the measurements are made only in dense regions, where galaxies may be over- or under- 
represented, depending on the cosmological model. Given the level of galaxy clustering in 
the SDSS, the plateau in M/L, for galaxies brighter than M r = —20, should be similar to 
the universal mean if <t 8 = 0.9, above the universal mean if a 8 > 0.9, and below the universal 
mean if a 8 < 0.9. As we will demonstrate in the next section, the value of a 8 for which 
the plateau is equal to the universal value depends on the bias of the galaxy sample under 
consideration, so for fainter luminosity thresholds the "unbiased" value of cr 8 is lower. 
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3.2. M/L of Clusters 

The M/L ratios shown in Figure 5, averaged over top-hat spheres of constant S, are not 
directly observable. In Figure 6, we make quantitative predictions for cluster sized halos that 
can be compared to observational data. For each a 8 , equation (4) is used to calculate M/L 
as a function of halo mass. Note that this calculation does not use simulations, so it can 
include galaxies down to M r = —18. We present predictions for galaxies with M r < — 18 and 
for galaxies with M r < —20, and refer to these mass-to-light ratios as M/Li 8 and M/L 2 o, 
respectively. 

Figure 6a plots M/L 18 , calculated using the soft central cutoff, against halo mass for all 
values of a 8 . Note that halo masses are proportional to Q m for fixed a 8 , so we list masses 
in wo.3^ _1 M and mass-to-light ratios in u , 3 hM Q /L Q . The M/L vs. M h curves derived 
from our modeling of w p (r p ) are similar in form to the parameterized function used by Yang 
et al. (2002) in their CLF modeling, and to the results of semi-analytic modeling by Benson 
et al. (2000). All the curves in panel (a) show a clear minimum at ~ 4 x 10 11 o;o.3 M Q . 
Below this minimum halo mass, M/L i8 increases rapidly, as it must to match the observed 
galaxy luminosity function to the steeper, low mass end of the halo mass function (see, e.g., 
Yang et al. 2002). At higher halo masses, M/L rises less rapidly and eventually reaches a 
maximum between 1 — 5 x 10 14 Uo^H^Mq, with higher values of a s reaching maximum 
M/L 18 at higher masses. At still higher masses, M/L 18 gradually declines, with the results 
for <j 8 = 0.6 falling by ~ 20% between 10 14 h' 1 M© and 10 15 h~ l M . 

The mean mass-to-light ratio of the universe, for galaxies brighter than M r < —18, can 
be calculated by integrating the Blanton et al. (2003) Schechter function fit to the observed 
r-band luminosity function of SDSS galaxies. This gives a luminosity density of 1.63 x 10 8 h 
L Q Mpc~ 3 , using an r-band absolute magnitude of the sun of 4.76, also taken from Blanton 
et al. (2003). Dividing the matter density, p m = Q m x p crit = 2.78 x 10 11 h 2 Vt m M & Mpc~ 3 , 
by this luminosity density gives the mean r-band mass-to-light ratio of the universe, 509 

/^Mq/Lq. 3 A similar value of 521 o;o.3 /iM /L is obtained by simply taking the number 
density in each magnitude bin listed in Zehavi et al. (2004b) Table 2 and using a discrete 
sum over all nine bins. The results in Figure 6a are similar to those of Figure 5a, but with 
a different model being unbiased with respect to (M/L); the halos above 10 14 hr l M for 
the a 8 = 0.8 model have M/L 18 & 490 lu .3 hM Q /L Q , which is very close to the cosmic mean 
value. For a 8 above and below 0.8, the cluster M/L 18 ratios are above and below this mean. 



3 To be precise, note that the Blanton et al. (2003) luminosity function is calculated in 01 r, the SDSS 
r-band redshifted to z — 0.1, and that this shifted bandpass is also used to define the luminosity threshold 
samples in Zehavi et al. (2004b). 
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Also plotted in Figure 6a are data from the Carlberg et al. (1996) analysis of clusters 
from the CNOC survey, with M/L in Gunn r-band calculated using mass estimates from 
the virial theorem. To properly compare our calculations to their results we have taken the 
values listed in their Table 4, which have been extrapolated to include all luminosities below 
their magnitude limit of —18.5 (a ~ 15% correction), and removed most of this correction 
up to our magnitude limit of —18.0 using the Schechter function parameters listed in their 
paper. This correction increases Mj L ratios by 6% relative to their stated values. We have 
also included a moderate correction for redshift evolution, since the mean redshift of the 
CNOC data is z = 0.3, while the SDSS data are centered on z = 0.1. Using the approximate 
correction factor of l0 015Az (Carlberg et al. 1996) we applied an additional 7% reduction 
of the luminosities. We ignore the slight differences between Gunn r and SDSS 01 r, which 
should be largely removed by the correction to solar units. The error-weighted mean of 
the Carlberg et al. data is 359 ± 32 /iM /L . 4 Inserting this value into equation (1), i.e., 
assuming that it represents (M/L), gives Q m ~ 0.2. However, Figure 6a demonstrates that 
the CNOC results are consistent with Q m = 0.3 if a 8 0.65, since the cluster M/L is then 
below the universal value. 

A number of authors have reported a trend of increasing M/L with cluster mass (e.g., 
Bahcall & Comerford 2002; Lin et al. 2004; Popesso et al. 2004), a result seemingly in 
conflict with the claimed plateau of M/L at high mass and with our results here. In the 
observational data and in our M/ L curves, there is a significant increase in M/ L from group 
masses (~ 10 13 /i _1 M ) to the cluster mass regime. At M > 10 14 h~ l M , the observational 
data of Bahcall & Comerford (2002) and Popesso et al. (2004) are consistent with a horizontal 
line. Lin et al. (2004) find a positive slope well into the cluster mass regime, but Kochanek 
et al. (2003), analyzing a similar sample with different methods, find a decrease in M/L 
with cluster mass, suggesting a significant systematic uncertainty in the detailed behavior at 
high masses. Our predicted M/L curves imply that a single relation between M and M/L 
is a poor approximation for any samples that extend below ~2x 10 14 h^ 1 M . Because of 
the observational uncertainties in the M/ L trend and the dependence of the predicted trend 
on the assumed value of a (see Figure 4a), we do not use this trend to draw cosmological 
conclusions. Instead, we use only the mean M/ L in the cluster mass regime, a quantity that 
is more robust observationally and theoretically (see Figure 4b). 

Panel (b) of Figure 6 shows the same calculation as panel (a), but now we only consider 
galaxies brighter than M r = —20. Using the Blanton et al. (2003) luminosity function, the 
mean M/L of the universe for this magnitude threshold is 923 cj .3 hM & /L Q . The predicted 



4 This mean value does not include the two clusters from the sample that show strong binarity. The error 
is computed by the bootstrap method. 
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M/L ratios still separate at high mass in the same proportions as in panel (a), but now the 
<7 8 for which (M/L) cl !=s (M/L) is 0.9 instead of 0.8 for M/L 1S , as with the numerical results 
presented in Figure 5. 

Figure 6c plots M/Li 8 against galaxy multiplicity rather than halo mass. We compute 
this relation by integrating over the halo mass function to calculate the contribution of each 
halo mass to the abundance of clusters at a given N gal , i.e. 

MIT IN i WlgVM m 

where all values with a subscript M are values at a given halos mass and dn/dM is the Jenkins 
et al. (2001) halo mass function. The probability P{N g&l \M) = (Akt>M gal ~V^at>Af/(jV gal _ i) 
for M > M min , since the distribution of satellite numbers in Poisson, and a halo with a satel- 
lite has a central galaxy by definition. Because of the exponential cutoff at the high-mass 
end of the halo mass function, there is asymmetric scatter into a given value of iV ga i; Poisson 
fluctuations around {N)m cause more halos of lower mass to scatter to high N ga \ than vice 
versa. This effect flattens out the curves relative to those of panel (a), but the asymptotic 
behavior at high multiplicity is similar. 

Figure 6d shows the same calculation for M/L 20 rather than MjL\%. At the same 
halo mass, the number of galaxies with M r < —18 is roughly five times the number with 
M r < —20, which is why the M/L curves reach their asymptotic values at N ga \ ~ 20 rather 
than ~ 100. As in panel (b), the asymptotic value of M/L for cr 8 = 0.9 is closest to the 
mean value. 

Figure 6d can be compared to the results from Bahcall et al. (2003a) analysis of clusters 
in the SDSS Early Data Release. We consider the cluster sample identified by the maxBCG 
method (see, e.g. Hansen et al. 2004), which characterizes richness by N-qcg, the number of 
galaxies close to the brightest galaxy in a restricted region of color-magnitude space. Bahcall 
et al. (2003a) report scaling relations 

LS. 6 (1O 1o L ) = 1.42^bcg, (6) 
o v (km s- 1 ) = 93AS, (7) 



for luminosity and velocity dispersion as a function of Nbcg- We have used the Schechter 
function parameters given in Bahcall et al. (2003a) to correct equation (6) from their observed 
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limit of M r = —19.8 to our threshold of M r = —20. The subscript 0.6 indicates that the 
luminosity is the total value within 0.6 h~ l Mpc of the cluster center, the radius at which all 
the cluster attributes are calculated by Bahcall et al. (2003a). The velocity dispersion can 
be converted to mass through a relation calibrated on gravitational lensing measurements of 
cluster masses (Bahcall et al. 2003b): 

M h (r < 0.6 h^Mpc) = 3.28 x lOV^ 67 ^ /i _1 M , (8) 

where ks is a small correction factor that depends on the mean overdensity of the halo within 
the defined radius (see Evrard, Metzler, & Navarro 1996). Since ks requires knowledge of 
the halo mass, we solve equation (8) by iteration, then combine with equation (6) to obtain 
M/L 2 o for 8 < A^bcg < 40, the range over which the scaling functions are valid. We 
convert from iV" B cG to the number of galaxies above M r = —20 using N gal = 0.14 A^q, an 
approximate scaling determined from the dependence of the mean cluster luminosity function 
on iV BC G (Wechsler et al. 2004). 

The shaded region in Figure 6d encloses the mean relation derived from equations 
(6) — (8) and is bounded by the 2-sigma errors on the scaling coefficient in equation (6). 
To make this comparison, we have also assumed that M/L 2 o does not vary from 0.6 hr 1 Mpc 
to the edge of the cluster. In our simulations, we find a modest increase in M/L from 0.6 
to 1.5 h~ x Mpc of about 20%, due to the bright central galaxy in each halo. The large 
contribution of the central galaxy to the overall luminosity of the cluster is also seen in the 
cluster luminosity functions of Hansen et al. (2004). However, the trend with radius is much 
smaller than the statistical errors on the scaling relations themselves. 

There are significant systematic uncertainties in our comparison because we combine 
scaling relations that have large individual uncertainties and intrinsic scatter. Future M/L 
measurements for larger SDSS cluster samples will enable more direct comparisons. For the 
current data, the mean relation plotted in Figure 6b seems consistent with the Carlberg et al. 
(1996) results. For a universe with Vl m = 0.3, the mean relation is consistent with a 8 = 0.6, 
while a lower value of tt m ks 0.12 is required to match the observations for a 8 = 0.9. 

At high halo masses, the values of M/L± 8 and M/L 20 are close to the universal values for 
<r 8 = 0.8 and a 8 = 0.9, respectively. The difference reflects the higher amplitude correlation 
function of more luminous galaxies. For a given value of a 8 , reproducing this trend requires 
putting a larger fraction of the more luminous galaxies in the strongly clustered, high mass 
halos. On large scales, the galaxy correlation function is ^ gg = b^ mm , where the galaxy bias 
factor b g is a number- weighted average of the halo bias factor bh(M): 
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where n g is the mean galaxy number density. Figure 7 plots b g for luminosity thresholds M r = 
— 18, —19, and —20, against the ratio (M/L) cl -=- (M/L), where (M/L) c \ is evaluated at M — 
5 x 10 14 h~ l M and (M/L) is the universal mass-to-light ratio. For a given luminosity 
threshold, b g and (M/L) c \ are decreasing and increasing functions of as, respectively, since 
matching the observed galaxy correlation function requires a lower bias factor for higher as, 
and a lower bias factor implies a smaller fraction of galaxies in high mass halos. At fixed a 8 , 
b g must be larger for high luminosity thresholds, and (M/L) c \ + (M/L) is correspondingly 
lower. In practice, we find in Figure 7 that (M/L) c \ ~ (M/L) for the value of a 8 that has 
b g ~ 1. There is no reason this must be exactly true, but our results are well captured by a 
simple rule of thumb: the cluster mass-to-light ratios for a given luminosity threshold and 
as are below the universal value if the large scale galaxy correlation function is positively 
biased, above the universal value if the correlation function is anti-biased, and equal to the 
universal value if the correlation function is unbiased. 



Attempts to measure M/L of cosmic structure often focus on galaxy clusters, since 
their masses can be estimated by the virial theorem, by more general dynamical models 
(e.g., Carlberg et al. 1997), by modeling their X-ray emission, or by weak gravitational 
lensing. Outside the virial region, there is still matter that is gravitationally bound to the 
cluster, but it is not in dynamical equilibrium, so the above methods (with the exception of 
weak lensing) are inapplicable. At the boundary of the infall region, where peculiar velocity 
cancels Hubble flow, the galaxy phase space density becomes infinite, creating caustic-like 
features in redshift space. The amplitude of these caustics is a measure of the escape velocity 
of the system (Kaiser 1987; Diaferio & Geller 1997; Diaferio 1999). The goal of the CAIRNS 
survey (Rines et al. 2003, 2004) was to identify these features and thereby measure the mass 
profile of clusters out to r ~ 10 hr x Mpc. 

To compare our predictions to the results of the CAIRNS survey, we identify cluster- 
mass halos (M > 3 x 10 14 h~ l M ) in our N-body simulations and calculate the average 
M/L ratio as a function of radius. Since the CAIRNS sample is taken from 2MASS K- 
band data, the M/L ratios need to be converted to r-band. The total luminosities listed in 
Table 8 of Rines et al. (2004) are in solar .fT-band units, corrected for incompleteness. Using 
r — K colors of 1.43 for the sun and 2.81 for elliptical galaxies (Pahre 1999), we multiply 



3.3. M/L in cluster infall regions 



-19- 



{M/L) K by io°- 4 ( 2 - 81 - L43 ) = 3.57 to get (M/L) r . We have multiplied each M/L ratio by 1.7 
to remove the luminosity from galaxies fainter than Mk = —22.81 for a proper comparison 
to our M r < —20 predictions. 

The CAIRNS M/L ratios have a specific geometry: the total light is in cylinders and 
the mass is in spheres. We calculate Mj L 20 ratios from our N-body simulations in the same 
way: centering on the most-bound particle of each halo with mass larger than 3 x 10 14 
h^ 1 M , the light at each projected radius r p is calculated within a cylinder that is extended 
10 h~ l Mpc in either direction from the cluster center, using the z-axis of the box as the 
line of sight. The total mass in dark matter particles is calculated in spheres of radius r p . 
The results for all a 8 values are shown in Figure 8. The ratio of these different geometries 
lowers M/L 20 , since the volume of the cylinder is larger than that of the sphere. In this 
spheres-on-cylinders calculation, the M/L 2 o values at 3 < r p < 10 hT 1 Mpc are lower than 
the mean, in contrast to the previous results. In tests that use equal cylinders for both light 
and mass, the Mj L 20 values, relative to the mean of the box, appear much as they did in 
Figure 6, with the a 8 = 0.9 model lying closest to the universal M/L at large r p . 

Points with error bars represent the CAIRNS measurements at r 200 (stars) and cluster 
turnaround radii (squares), taken from Table 8 of Rines et al. (2004) and converted to r-band 
as described above. Here r 2 oo is the radius at which the overdensity of the cluster is 200 
times the critical density (666 times the mean density assuming Q m = 0.3), which is close 
to 1 h^ 1 Mpc for all the clusters in the CAIRNS sample. The error-weighted mean of the 
r 2 oo results is 424 ± 37 hM Q /L & . For Q m = 0.3, this value for (M/L 20 ) c i implies a 8 0.7, 
similar to the CNOC results and the SDSS results shown in Figure 6. For cr 8 0.9, a 
value of Q m ~ 0.18 is required to bring the numerical results into agreement with the data. 
The error- weighted mean of the Mj L 20 values from the cluster -l-infall regions is significantly 
lower, 313 ± 32 /iM /L . For a 8 = 0.9, this ratio implies Vt m m 0.14. 

The discrepant conclusions between the virial and infall mass estimates suggest that 
the mass profiles inferred by Rines et al. (2004) are steeper than those predicted by our 
simulations. To investigate this point further, we plot in Figure 9 the ratio of the mass 
M tot within the turnaround radius r max to the mass M 2 oo within r 200 , as a function of r max . 
The observationally inferred mass ratios lie in the range 1.2 to 2.2, with a trend of larger 
mass ratios for increasing r max . Solid curves show the results for our simulated clusters. 
While M tot and M 2 oo in our simulations both scale with fi m , there is still a small dependence 
on Q m because we choose clusters above 3 x 10 14 h~ 1 M Q , and therefore select a different 
sample for different Q m . Note that the observational data points are independent of the 
galaxy luminosity measurements and the theoretical curves are independent of our galaxy 
bias models, since there ratios refer to mass alone. The theoretical curves lie above the data 
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points by ~ 50%, similar to the difference between the average virial and infall mass-to- light 
ratios. Our best guess is that the caustic method systematically underestimates the infall 
masses by ~ 50%, but it is of course possible that the observationally inferred mass profiles 
are correct and conflict with the generic predictions of the ACDM model. 



We have examined mass-to-light ratios of large scale structure in cosmological models 
that are constrained to match observed real space galaxy clustering. Specifically, we consider 
models in which the shape of the linear matter power spectrum is held fixed and the galaxy 
halo occupation distribution is adjusted to reproduce Zehavi et al.'s (2004b) measurements of 
the projected galaxy correlation function. For power spectrum normalization <7g in the range 
0.6 — 0.95, we are able to find HOD parameters that yield acceptable fits to the observed 
w p (r p ), even with a restricted, 3-parameter HOD prescription. For each value of <r 8 , the M/L 
ratios in high overdensity regions are approximately independent of top-hat smoothing scale 
in the range 1.5 — lOh^ 1 Mpc, and the M/L ratios of virialized halos climb from a minimum 
at Mh ~ several x 1O 11 /iM /L q to an approximately flat plateau in the cluster mass regime. 
However, this plateau only corresponds to the true universal Mj L for a particular choice of 
<7g, the one for which the large scale galaxy correlation function is unbiased. One therefore 
cannot take the existence of a plateau in M/L as a function of scale or of halo mass as 
evidence that one has measured the universal M/L. Estimates of Q m that multiply cluster 
mass-to-light ratios by the observed luminosity density make the implicit assumption that 
the galaxy distribution is unbiased. 

Given the SDSS clustering measurements, we expect cluster mass-to-light ratios to be 
representative of the universal value for cr 8 pa 0.8 if one is considering galaxies with M r < —18, 
or for (j 8 pa 0.9 if one is considering galaxies with M r < —20. For lower a 8 , cluster M/L ratios 
lie below the universal value because galaxies must be overrepresented in dense regions to 
match the observed clustering. Conversely, cluster M/ L ratios lie above the universal value 
for higher ct 8 - 

Averaging our results over the same masses as the CNOC cluster sample, our results 
are well described by the relations 



4. Implications and Outlook 




(10) 
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for luminosity thresholds of M r = —18 and M r = —20, respectively. As discussed in §2, we 
estimate that uncertainties in the HOD fits introduce a ~ 10% systematic uncertainty in 
the normalization of these relations, though our present investigation is not exhaustive. The 
luminosity function of SDSS galaxies implies universal mass-to-light ratios of (M/L 18 ) = 
509cjo.3 /iM /L and (M/L 2 q) = 923cjo.3 ^M /L , where tuo.3 = fi m /0.3. The statistical 
uncertainty in Blanton et al.'s (2003) luminosity density estimate implies a 2% uncertainty 



If we combine equation (10) with the mean mass-to-light ratio of CNOC clusters, 
(M/Li 8 ) c i= 359 ± 32/iM /L , we obtain the constraint 



or o"8^m 6 = 0.33 ± 0.03. Here we have added in quadrature the 9% statistical error on the 
mean cluster M/L, the 2% statistical error in the mean luminosity density, and our estimated 
10% systematic error in the normalization of equation (10), but we have not considered other 
possible sources of systematic error. We find similar results comparing to the Mj L ratios of 
SDSS clusters inferred from Bahcall et al. (2003a) or to the M/L ratios of the virial regions 
of clusters found by Rines et al. (2004), but the systematic uncertainties are larger and 
harder to quantify because of our reliance on mean scaling relations in the former case and 
the complications of passband and geometry conversions in the latter. If we used the Rines 
et al. (2004) M/L estimates from cluster infall regions, we would infer a somewhat lower 
value of cr 8 f^ 6 . 

Our estimate (12) derived from cluster mass-to-light ratios conflicts with recent es- 
timates obtained by combining CMB anisotropy measurements with the large scale galaxy 
power spectrum, Type la supernova data, the Lya forest flux power spectrum, and weak grav- 
itational lensing, which tend to favor a 8 0.9, Q m ~ 0.3 (e.g., Spergel et al. 2003; Tegmark 
et al. 2004; Seljak et al. 2004). For example, Seljak et al. (2004) quote a 8 = 0.897:^03?, 
2811q;o2i as their combined constraint for a 6-parameter, spatially flat ACDM model, 
corresponding to <J 8 Vt^ = 0.419 ± 0.026. Our conclusion agrees well with that of van den 
Bosch et al. (2003), who find that "concordance" values of <j 8 = 0.9, Vt m = 0.3 are favored 
by their conditional luminosity function analyses of the 2dFGRS only for cluster mass-to- 
light ratios of 75O/iM /L . If they impose a more observationally plausible constraint of 
(M/L) cl = 350 ± 7O/iM /L (in 6j-band), they find Vt m = 0.25l£$ and a 8 = 0.78 ± 0.12, in 
good agreement with equation (12). Although our calculation is similar in overall concept 
to van den Bosch et al.'s, we parameterize the problem in a completely different way, fit 
different constraints in a different order, use different approximations, and analyze measure- 



in (M/L). 




(12) 
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ments from an independent galaxy redshift survey, red-selected instead of blue-selected. In 
contrast to our fixed linear P(k), van den Bosch et al. (2003) vary the power spectrum shape 
parameter linearly with Q m . The agreement of the two results is therefore a good indication 
that the conclusions are robust to details of the measurements or analysis procedures. 

Equation (12) is nearly identical to the constraint 0" 8 f2^ 6 = 0.33 ±0.03 obtained by Bah- 
call et al. (2003b) in their analysis of the mass function of clusters in the SDSS Early Data 
Release. Earlier analyses of cluster mass functions have generally yielded higher normaliza- 
tions of this constraint (e.g., White et al. 1993; Eke et al. 1996; Bahcall & Fan 1998; Henry 
2000). While the mass function and mass-to-light ratio methods both incorporate cluster 
masses, they are physically distinct: the former has no dependence on galaxy luminosities, 
and the latter uses an average mass and is therefore insensitive to scatter in observational es- 
timates. Low values of a 8 or Q m would also help explain observational estimates of the galaxy 
pairwise velocity distribution, which appear to conflict with predicted values for <r 8 = 0.9, 
Q m = 0.3 (Yang et al. 2004). 

Our HOD modeling assumes that satellite galaxies in halos have a radial profile cor- 
responding to an NFW model with the concentration predicted for VL m = 0.3. To test our 
sensitivity to this assumption, we lower all of the assumed halo concentrations in the as = 0.9 
model by 30% and refit w p (r p ), which results in slightly different P(N\M) and thus slightly 
different M/L. We find a 1.8% difference in (M/L) c \ for this HOD model, so radial pro- 
files are not an important source of systematic uncertainty. As noted in §2, changing the 
assumed linear matter power spectrum from the Efstathiou et al. (1992) parameterization 
to a CMBFAST calculation makes a significant difference to the x 2 values of w p (r p ) fits 
but minimal change to the HOD parameters themselves, so our M/L predictions are not 
sensitive to modest changes in the power spectrum shape. The main source of systematic 
uncertainty is our (M/L) cl predictions is therefore the 10% uncertainty associated with the 
HOD parameterization and the w p (r p ) measurements, discussed at the end of §2 (see Fig. 
4). 

Turning to the observational uncertainties, cluster Mj L ratios could be underestimated 
if masses are biased low or luminosities are biased high. Cluster mass estimation is a chal- 
lenging problem, but the generally good agreement between virial mass estimates, Jeans 
equation estimates, and estimates from X-ray or weak lensing data (see Bahcall & Comer- 
ford [2002] and references therein) argues against a systematic error as large as a factor of 
588/359 = 1.64, the ratio of our predicted M/L for Q m = 0.3, a 8 = 0.9 to Carlberg et al.'s 
(1996) observational estimate. Luminosities could be biased high if background subtraction 
methods do not adequately account for overdense structures surrounding clusters, an issue 
that warrants further investigation with realistic mock galaxy catalogs. We see no obvious 
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holes in the CMB + large scale structure analyses, but in contrast to the approach taken 
here and in van den Bosch et al. (2003), these methods of inferring cosmological parameters 
rely on a detailed theoretical model of primordial fluctuations and their linear evolution, 
and they are sensitive to quantities like the electron scattering optical depth, the CMB 
tensor-to-scalar ratio, and the curvature of the inflationary fluctuation spectrum. 

At present, it seems plausible that the results of these different methods can be rec- 
onciled without major revisions, if og and fl m lie at the low end of the ranges allowed by 
Tegmark et al. (2004) or Seljak et al. (2004), cluster M/L ratios are somewhat higher than 
the estimates shown in Figures 6 and 8, and our HOD parameterization leads to a modest 
over-prediction of M/L for a given as and Q m . Fortunately, the remaining uncertainties 
can be substantially reduced in the near future. The SDSS redshift survey is now large 
enough that rich clusters can be identified directly from the redshift survey itself, reducing 
(though not eliminating) problems of contamination and background subtraction. Follow-up 
observations of these systems can provide consistency checks of mass estimates via galaxy dy- 
namics, X-ray modeling, and weak lensing. The HOD parameter constraints can be greatly 
improved by bringing in additional clustering measurements, most notably the group multi- 
plicity function (A. Berlind et al., in preparation) and the projected three-point correlation 
function. Redshift-space distortion analysis can yield independent constraints on cr 8 and Vt m . 
Advances on all of these fronts should soon show whether the current tension in parameter 
estimates arises from an accumulation of systematic errors or instead signals the need for a 
new physical ingredient in the standard cosmological scenario. 
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A. Halo and Galaxy Bias 



Halo clustering is biased relative to that of the underlying mass distribution, by an 
amount that depends on halo mass. Halo bias factors are an important ingredient in analytic 
calculations of galaxy clustering, including those in this paper. Following the pioneering 
work of Mo & White (1996), a number of authors have investigated halo bias using N- 
body simulations (e.g., Porciani et al. 1999; Sheth & Lemson 1999; Sheth & Tormen 1999; 
Jing 1998, 1999; Sheth, Mo, & Tormen 2001). However, many of these studies are based 
on simulations of either low mass resolution or limited box size (the main exception is the 
recent study of Seljak & Warren [2004], which we discuss below). Most previous studies 
have also compared different cosmologies such as open and standard CDM to one version 
of ACDM, rather than focusing on variants of the ACDM cosmology with a wide range of 
o" 8 as we have done here. We have also performed an identical set of simulations with a 
power spectrum shape parameter of T = 0.12, which significantly increases the large-scale 
power and reduces the power on small scales relative to our standard choice of T = 0.2. 
Our simulations are therefore well suited to investigate halo bias for the cosmological models 
of the greatest interest today, and to investigate the dependence of bias factors on power 
spectrum shape or normalization. The use of five realizations and a reasonably large volume 
(253 h^ 1 Mpc on a side) yields good statistics for high mass halos. Note that we define halos 
using a friends-of-friends algorithm with linking length I = 0.2n -1 / 3 . Alternative definitions 
would yield slightly different halo masses (see Hu & Kravtsov 2003) and would therefore 
require a slightly different formula for bh{M). 

The halo bias factor can be defined by the ratio of halo and mass autocorrelations, or 
power spectra, or using the halo-mass cross-correlation. Since we are interested in modeling 
galaxy autocorrelations, we adopt the definition b\(M) = ^(r, M)/£ m (r), where ^ h (r,M) is 
the autocorrelation function of halos of mass M and £ m (r) is the non-linear matter correlation 
function measured from the simulations. Sheth et al. (2001, hereafter SMT) give an analytic 
formula for b h , motivated by the analytic model of Sheth & Tormen (1999) but empirically 
calibrated against numerical simulations: 



b h {u) = 1 + 



^{av 2 ) + ^Jab(av 2 ) l - C - 



2\c 



W 2 ) c + 6(1 -c)(l -c/2) 



(Al) 



where 5 C = 1.686 is the critical overdensity required for collapse and v = 5 c /a(M), with 
a(M) the linear theory rms mass fluctuation in spheres of radius r = (3M/47rp) 1//3 . The 
three parameters in this equation are a = 0.707, b = 0.5, and c = 0.6, as listed in SMT. 



We divide the halos in our simulations into bins separated by factors of two in mass. The 
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logarithmic center of the lowest bin is 1.22 x 10 12 (f2 m /0.3) h 1 M , below M* for <j 8 > 0.6. 
We calculate halo bias by averaging 

^Jih/im for radii 4 < r < 12 hr x Mpc, a regime in which 
the ratio of the halo and matter correlation functions is approximately constant and noise is 
not a factor. The results are plotted in Figure 10. The dashed line, which shows the SMT 
bias relation, is significantly higher than the values of bh calculated from the simulations for 
both values of V. A better fit to the calculations is shown with the solid line, which plots a 
bias relation of the same form as equation (Al), but with a = 0.707, b = 0.35, and c = 0.80. 
With these parameter values, the formula gives accurate fits to our numerical results for 
the full range of a 8 values, and it works equally well for T = 0.2 and T = 0.12. Increasing 
either the lower or upper bounds of the radial range over which bh is calculated does not 
appreciably change the bias values or the quality of the fits. 

The inset box in Figure 10 plots an example of using equation (9) to calculate the galaxy 
bias with our modified parameters of equation (Al), which we compare to the galaxy bias 
in the simulations calculated by the same method as the halo bias for five values of a s . The 
plot symbols represent the galaxy bias calculated from the simulations, and the solid lines 
represent the analytic calculations (eq. [9]) using equation (Al) with our new parameters. 
For this test, the HOD parameters used to populate the simulations are taken from Table 
1 of Tinker et al. (2005), parameters similar to the M r < —20 sample. The analytically 
calculated bias factors differ from the simulation results by < 1%. 

Recently Seljak & Warren (2004, hereafter SW) proposed a new empirically determined 
halo bias relation, empirically calibrated on large simulations (up to 768 3 particles) with 
cosmological parameters close to the best-fit values from CMB and large-scale structure 
measurements. The dotted line in Figure 10 shows the SW formula evaluated for T = 0.2 
and <t 8 = 0.9, which fits our numerical data accurately. However, the SW formula is expressed 
in terms of M/M*, where cr(M*) = 5 C , instead of v — 5 c /a(M), and the mapping between 
M/M* and v depends on the amplitude and shape of the power spectrum. We find that the 
SW formula becomes a poor fit to our results at low as, and the discrepancies are worse for 
T = 0.12, as illustrated in the inset box by calculating the galaxy bias parameters with the 
SW halo bias formula. Thus, their formula is accurate for models close to the cosmological 
concordance model, but our modified version of the SMT formula, with v as the halo mass 
parameter, applies more universally. 

We have not investigated alternative definitions of bias using the power spectrum or 
mass cross-correlation, or using linear instead of non-linear matter clustering. We also note 
that our numerical results do not extend much below v — 1, and simulations of smaller 
volumes or larger dynamic range are needed to test equation (Al) in the low mass regime. 
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B. The Analytic Model 

Our analytic calculation of the two-point galaxy correlation function is similar to that 
presented by Zheng (2004) and used in Zehavi et al.'s (2004a,b) modeling of the SDSS 
w p (r p ). Here we report improvements in the procedure that have been incorporated in our 
present analysis. For completeness and clarity we describe the method from start to finish. 
The new ingredients are the use of the modified halo bias formula of Appendix A and the 
more accurate treatment of halo exclusion described by equations (BIO) and (B12). The 
Zheng (2004) procedure was calibrated and tested using the 144 h^ 1 Mpc, 256 3 particle GIF 
simulation of Jenkins et al. (1998), which was also used to calibrate the SMT formula. Our 
use of the five 253 hr x Mpc, 360 3 particle simulations described in §3 allows us to achieve a 
more accurate calibration and to test the procedure for a range of cosmological parameters. 

The correlation function is defined as the probability above random of there being a pair 
of objects at separation r. In the HOD context, a pair of galaxies can reside within a single 
halo or come from two distinct halos. These two contributions are computed separately and 
combined to get the full correlation function, i.e. 



£(r) = [l+U(r)]+6h(r). (Bl) 

The "1+" arises because it is the pair counts, proportional to 1 + £ih and 1 + £ 2 h) that sum 
to give the total pair counts, proportional to 1 + £. The one-halo term is calculated in real 
space through (Berlind & Weinberg 2002) 

^+f (r) 1 r iM dn ^ n{n ~ 1))m 1 f'( r \ mi) 

1 + Uir) -^J dM dM 2 2^(M) F feJ' (B2) 

where n g is the mean number density of galaxies, dn/dM is the halo mass function (Sheth & 
Tormen 1999; Jenkins et al. 2001), and (N(N — 1))m/2 is the average number of pairs in a 
halo of mass M. The function F(x) is the average fraction of galaxy pairs in a halo of mass 
M (or virial radius i? v ir) that have separation less than r, which is related to the halo density 
profile, p m (r), and F'(x) is its derivative. In practice, F(x) must be treated differently for 
central-satellite galaxy pairs and satellite-satellite pairs. In the former, the pair distribution 
is proportional to the volume- weighted density profile, F'(x) oc p m (r)r 2 , normalized to one. 
For the latter it is derived from the halo profile convolved with itself, a calculation that can 
be done analytically for an NFW profile (Sheth et al. 2001). The average number of one-halo 
pairs in the range (x, x + dx) in halos of mass M can be written explicitly as 
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F'(x) dx = {N cen N SSLt } M F^ s (x) dx + 
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K( x ) dx > 



(B3) 
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where the subscripts cs and ss refer to central-satellite pairs and satellite-satellite pairs 
respectively. For a Poisson distribution of satellite occupation, {N sat (N sat — 1)) = (iV sat ) 2 . 
We use an NFW profile with concentration parameters as a function of halo mass calculated 
by the method of Bullock et al. (2001) and Kuhlen et al. (2004). 

The one-halo term dominates £ 3 (r) at small scales, while the two-halo term fully accounts 
for all galaxy pairs at separations > 5 h^ 1 Mpc. The transition region between one-halo 
and two-halo dominance is difficult to model because only certain regions of the halo mass 
function can contribute to the two-halo term at small scales. The range of halo masses 
included must ensure that halo pairs do not overlap, since such halo pairs would be merged 
into a single halo by the friends-of-friends scheme that we use to define halos in the first place. 
It is this one-halo to two-halo transition that causes w p (r p ) to deviate from a power-law at 
scales near 1 fir 1 Mpc (Zehavi et al. 2004a,b); as r increases, £ih drops rapidly, while the 
rise in £ 2 h is regulated by halo exclusion. For brighter galaxies, which preferentially occupy 
high-mass halos, the rise in £ 2 h occurs at a larger r, making the deviation from a power-law 
greater for brighter galaxy samples. 

Since the radial distribution of galaxies within halos must be accounted for in the calcu- 
lation of the two-halo term, the calculation itself is done in Fourier space, where the convolu- 
tions with the halo density profile become multiplications instead (Scherrer & Bertschinger 
1991; Seljak 2000; Soccimarro et al. 2001). In the method of Zheng (2004), halo exclusion 
was treated by only including halos with virial radii less than half the value of r for which 
£ 9 (r) is being calculated, i.e. with masses below 



where A is the virial overdensity of the halo, relative to the mean density, which we have 
chosen to be 200. With this implementation of halo exclusion, the calculation of the two-halo 
term in Fourier space is 




(B4) 




-i 2 



P*\k,r) = P m (k) 



dM—(N) M b h (M,r)y g (k,M) 



(B5) 



where y g (k, M) is the Fourier transform of the halo density profile (e.g. Cooray & Sheth 
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2002), bh(M,r) is the halo bias at separation r, and the restricted number density n' g is the 
average number density of galaxies that reside in halos with M < M lim , 

f Mlim dn 

n' 9 = J o rfM-(JV> tf . (B6) 

The matter power spectrum, P m (k), is the non-linear form given by Smith et al. (2003). At 
large r, equation (B5) can be thought of as simply multiplying the non-linear matter power 
spectrum by the galaxy pair-weighted halo bias factor to obtain the galaxy power spectrum. 
At smaller separations the finite size of the halos must be taken into account. The scale 
dependence of halo bias also becomes important at r < 3 h^ 1 Mpc. Parameterizing the 
scale dependence by the amplitude of the non-linear matter correlation function, the scale 
dependence of halo bias is well described by 

b 2 (M r) = b 2 (M ) [l + l-17em(r)] L49 

1 ' ' K ] [l + 0.69£ m (r)] 2 - 09 ' 1 ' 

where b(M) is the large-scale bias, for which we have used the bias relation given in Appendix 
A. Equation (B7), determined from inspection of our numerical simulations, is fairly accurate 
for the full range of <t 8 values and both values of T explored in this paper. 

At a given r, equation (B5) is solved for all fc, then converted to real space by 

^ = ^ r p ^ k ' r)k2 ^ Ldk ' <B8) 

where £,' 2h (r) denotes that we have calculated the two-halo term for a restricted range of the 
halo mass function. The value calculated in equation (B8) is converted to a probability over 
random for the entire halo (and therefore galaxy) population by 

l + 6h(r)=(^ 2 [l + ekW]. (B9) 

The virtue of this implementation is that the integral over mass in equation (B5) is 
calculated once and squared, instead of being a double integral over different halo pairs. 
This approximation, however, neglects galaxy pairs from halos larger than M lim paired with 
smaller halos. In Figure 11, we compare this analytic method to numerical results. For 
this comparison we use the correlation functions calculated from the N-body simulations 
described in §3, but following Tinker et al.'s (2004) practice of drawing satellite galaxy 
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populations from the appropriate NFW profile instead of randomly sampling the friends- 
of-friends dark matter distribution, as done in §3. We use the ag = 0.8 output with HOD 
parameters M min = 1.11 x 10 12 , M x = 2.53 x 10 13 , and a = 1.01. 

Figure 11a compares the N-body two-halo term to equation (B9) calculated both with 
our new halo bias parameters and with the original parameters of the SMT function. At 
r > 10 h^ 1 Mpc, the original SMT relation over-predicts the galaxy bias by ~ 15% and the 
correlation function by ~ 30%. The new bias function yields an excellent match at these 
scales. At smaller scales, this method of halo exclusion under-predicts the number of two- 
halo pairs, with errors greater than 50% at 1 h~ l Mpc. Although the one-halo term begins to 
dominate at this scale, the large error in £211 is still apparent in the total correlation function, 
shown in panel (b). The ~ 5% error at the smallest scales is due to the fact that the halo 
mass function of our simulations is not precisely represented by the Jenkins et al. (2001) 
function assumed in the analytic calculation. The two-halo term is much less sensitive to 
the mass function, and is not affected by this difference. 

The numerical test in the Appendix of Zehavi et al. (2004a) showed a much smaller 
discrepancy on large scales because the 144 h' 1 Mpc GIF simulation used for the test (and 
for the calibration of the SMT bias factors) has a high amplitude of large scale clustering 
for low mass halos. With our multiple, larger volume simulations, the need for lower halo 
bias factors is evident, and these in turn drive the need for a more accurate treatment of 
halo exclusion. Zehavi et al. (2004a,b) find a low x 2 /d.o.f. fitting the projected correlation 
function of luminous, M r < —21 galaxies, while we find a relatively high x 2 /d.o.f. for this 
sample (see Table 1). However, as noted in §2, the x 2 values of w p (r p ) fits are sensitive to the 
difference between a T = 0.2 power spectrum and a (presumably more realistic) CMBFAST 
power spectrum, even though the best-fit HOD parameters are not. Combination of our 
present £ 9 (r) calculation with the CMBFAST power spectrum produces a similar w p (r p ) 
for M r < —21 galaxies, with similarly low X \ to that obtained with the Zheng (2004) 
prescription and a T = 0.2 power spectrum. Thus, this more accurate modeling leads to the 
same bottom line conclusion as Zehavi et al. (2004a,b), and similar HOD parameters. 

We now return to the halo exclusion problem. Under the assumption of spherical halos, 
all two-halo pairs would be accounted for by summing all the galaxies from halo pairs for 
which the sum of virial radii is smaller than the separation, i.e. -R v iri + -R v ir2 < t. For this 
"spherical halo" exclusion, equation (B5) must be modified to 
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J dM 2 — {N) M2 b h (M2,r)y g (k,M 2 ), 



/•M lim , 2 dn 



(BIO) 



where Mi im l is the maximum halo mass such that i? v i r (Mi im l ) = r — i? vir (M min ) and Mi im 2 
is related to Mi by R V u{M\ iui ,2) — r ~ Rvh(Mi). Since the upper limit of the second in- 
tegral depends on the integrand of the first, equation (BIO) must be solved as a double 
integral, making it more computationally expensive than equation (B5), but significantly 
more accurate by increasing the number of small-separation two-halo pairs counted. 

The range of halo masses over which the two-halo term is calculated is different than 
that in equation (B5), and the restricted number density n' g must reflect that change. Using 
this new method of halo exclusion, n' g becomes 



Figures 11c and lid compare our numerical results to the spherical exclusion method. At 
r = 1 h~ x Mpc, where the previous method resulted in a ~ 50% error, the error has been 
reduced to ~ 25%. At r = 2 h' 1 Mpc, the ~ 20% error in the previous method is eliminated 
entirely. 

Halos are not spherical objects, however. They are triaxial objects that can exhibit 
significant flattening (see, e.g. Jing & Suto 2002). This can lead to halo pairs which are 
closer than the sum of their virial radii. By assuming a lognormal distribution of ellipticities 
with mean flattenings motivated by simulation results, one can determine the probability that 
halos of a given separation are allowed. This probability, P(x), where x = r/(i? virl + it! vir2 ) , 
can be used to modify the two-halo calculation to 



where the limits in equation (B12) are both infinity (but in practice can be cut off at some 
reasonably large value) but the calculation is still a double integral because the value of Mi 
is used in the ellipsoidal exclusion probability in the integral over M 2 . 

We investigated the ellipsoidal exclusion probability by a Monte Carlo approach. We 
assume a lognormal distribution of halo axis ratios q with dispersion 0.2 and means = 




(BH) 




(B12) 
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b/a — 0.9 and q c = c/a = 0.8, motivated by the results of Jing & Suto (2002), then assume the 
axis ratios have a lognormal distribution in q with a dispersion of 0.2. By randomly selecting 
ellipticities and orientation angles for halos of a given mass ratio and separation, we find 
that the probability of non-overlapping halos is well approximated by P(y) = (3y 2 — 2y 3 ) 
in the range < y < 1, where y = (x — 0.8)/0.29. At y < 0, P{y) = 0, and at y > 1, 
P(y) = 1. The restricted number density, n' g , must also be calculated in this way. Equation 
(B6) becomes 



The results of the ellipsoidal halo exclusion are presented in panels (c) and (d) of Figure 
11. The ellipsoidal exclusion approach is an improvement over spherical exclusion; at r = 1 
hr l Mpc, the error in the two-halo term has been reduced to ~ 10%. Choosing more extreme 
values of qb and q c does not significantly change the results. This ellipsoidal exclusion method 
is the one we have used for the calculations in this paper. 

Although this approach is more accurate, the added computation time can become 
prohibitive when fitting observed data to high precision, where many hundreds of iterations 
are required. We have therefore created a halo exclusion approach that mimics the results of 
ellipsoidal exclusion but can be written as separable integrals in the Fourier-space calculation 
of P|g (k, r). Although the evaluation of equation (B12) is CPU intensive, the calculation of 
the restricted number density under the ellipsoidal approach is relatively rapid. For a more 
efficient scheme, we recalculate the ellipsoidal n' g first, and then Mu m in equation (B6) is 
increased until the restricted number density matches that of the ellipsoidal calculation. We 
then use equation (B5) to calculate £211- 

The results of this approximation, which we call n^-matched, are shown in panels (c) 
and (d) of Figure 11. The loss of accuracy relative to the full ellipsoidal treatment is minimal. 
In further tests with multiple values of as and HOD parameters that give higher and lower 
galaxy space densities, we find similar results. 

Although the methods introduced here significantly improve the analytic calculation of 
£, g {r) for specified cosmological and HOD parameters, there is room for further investigation 
and improvement. Outstanding issues include the small but non-negligible dependence of 
the halo mass function shape on cosmology (i.e., the non-universality of the Jenkins et al. 
[2004] formula), the effect of scatter in concentrations and halo ellipticity on the one-halo 
term, bias factors of low mass halos, effects of cosmology on the scale-dependence of halo 
bias, the dependence of the approximation on the halo definition, and the interaction of 
all of these effects with the treatment of halo exclusion. We are presently investigating a 
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number of these issues. The long-term goal is to ensure that errors in the calculation of 
£ g (r) for specified parameters are a negligible source of uncertainty in the inference of HOD 
and cosmological parameters from observational data. Because of the high precision of the 
clustering measurements, it is not clear that we have yet reached this goal. 
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Fig. 1. — (a) The HOD fits are compared to the w p (r p ) data for galaxy samples M r < —20 
(lower points) and M r < —21.5 (upper points). The five solid curves represent the five values 
of cr 8 , with the lowest curve being cr 8 = 0.6 in both fits, and each subsequent curve going 
in order of increasing cr 8 . (b) The projected dark matter correlation functions for Q m = 0.3 
and all five values of erg- The curves, from lowest to highest, go in order of increasing <7 8 . 
For comparison, the w p (r p ) data are plotted as well. 
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Fig. 2. — The conditional luminosity function, <&(L\M h ), is plotted for three different halo 
masses, log (M h/M Q ) = 12.5, 13. 5, a dn 14.5, and for all five values of <r 8 . The curves are all 
normalized such that the area under each curve is (N)m for M r < —18. The dashed line in 
each plot is an example of the contribution to $(L|M^) by satellite galaxies only (here for 
a 8 = 0.95). At M = 3 x 10 12 h' 1 M , satellites contribute only to the faint end of $(L|M h ), 
while for the cluster-mass halos of panel (c), satellites dominate the luminosity function for 
all but the brightest magnitude bin. 
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Fig. 3. — A comparison between soft and hard central galaxy cutoffs. Panel (a) demonstrates 
that the different central occupation functions do not change the predicted w p (r p ). The open 
squares are the data for M r < —20 galaxies from Zehavi et al. (2004b). The solid line is 
the HOD fit to the data using a hard cutoff. The solid circles represent the same HOD 
parameters, but w p (r p ) is recalculated for a soft cutoff. Panel (b) shows the difference in 
$(L|M/j) from the two methods for M = 3 x 10 13 h~ l M . A hard cutoff, shown with the 
solid histogram, produces a significant bump at M r ~ —22, which is smoothed out by the 
soft cutoff, shown by the dotted line and filled circles. Panel (c) shows the M/L ratio as 
a function of halo mass for the two methods. The solid line represents the hard cutoff and 
the dotted lines represents the soft cutoff. The differences are small, and nearly negligible 
at high masses. 
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Fig. 4. — (a) M/L ratio as a function of halo mass for five values of a: 0.9, 0.95, 1.0, 1.05 and 
1.1. All calculations assume a s = 0.9. The shaded region is the full-width at 10% maximum 
of the distribution of halo Mj L inferred from a flexible HOD parameterization with no fixed 
value of a. See text for further details, (b) The change in the mean cluster M/L ratio as a 
function of a relative to a = 1. The mean is calculated for the same distribution of cluster 
masses as the CNOC cluster sample. 
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Fig. 5. — The mean M/L ratio in spheres of radius r as a function of the enclosed overdensity 
1 + 5, for r = 1.5, 5.0, and 10.0 h^ 1 Mpc (panels a, b, and c, respectively). Results are 
calculated from the numerical simulations and refer to luminosities of galaxies with M r < 
—20. Different symbols represent different a s values with the same coding as Figure 2; at 
high 5, the as = 0.6 curve is always lowest. 
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Fig. 6. — Mass-to-light ratio (r-band) as a function of halo mass (top panels) or richness 
(bottom panels). Luminosities are for galaxies brighter than M r = —18 (left) and M r = —20 
(right). From bottom to top, curves represent <7 8 = 0.6, 0.7, 0.8, 0.9, 0.95. Dashed horizontal 
lines represent the mean M/L of the universe. In panel (a), open squares with error bars 
are the CNOC data of Carlberg et al. (1996). In panel (d), the shaded region represents the 
M/L ratio of SDSS clusters based on Bahcall et al. (2003a). 
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Fig. 7. — Galaxy bias, calculated for the SDSS samples from our HOD models, as a function 
of mass-to-light ratio at M = 5 x 10 14 o;o.3 M Q , relative to (M/L). The different plot 
symbols represent three different galaxy magnitude thresholds, M r = —18,-19, and —20. 
For each sample, the data points represent the five different a 8 values used. The data always 
follow the monotonic trend of increasing a 8 for decreasing b g (e.g., <t 8 = 0.6 at the far left 
and <t 8 = 0.95 at the far right). 
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Fig. 8.— Mean M/L ratios of clusters (with M > 3 x 10 14 ftr 1 M ) as a function of 
projected separation from the cluster center, calculated from the numerical simulations, 
from bottom to top, the curves represent <t 8 = 0.6,0.7,0.8,0.9,0.95. Following Rines et 
al. (2004), we calculate masses in spheres and luminosities (for M r > —20) in truncated 
cylinders, which depresses M/L values below those shown in Figures 6b and 6d. Stars 
and open squares represent the M/L estimates of Rines et al. (2004) at r 2 oo and cluster 
turnaround radii, respectively, converted from i^-band to r-band as described in the text. 
The dashed horizontal line shows the mean M/L20 of the box. 
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Fig. 9. — Ratio of the total mass within r max to the mass within r 2 oo from the Rines et al. 
(2004) data (open circles) and our simulations (curves). The predictions depend slightly on 
the value of Q m assumed in calculating cluster masses, since we consider only halos with 
M > 3 x 10 14 /j _1 M Q . The curves represent Q m = 0.1,0.3, and 0.45, and are calculated for 
o"s = 0.8; results for other as are similar. 
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Fig. 10. — The halo bias factor is shown as a function of v for the five different values of 
erg for both T = 0.2 and T = 0.12. The values of the bias factor were calculated from the 
asymptotic value of the correlation function at large scales: = >j£h/£m- The T = 0.12 
results are offset by 0.2 dex to avoid crowding. The dashed line is the bias relation given by 
Sheth, Mo, & Tormen (2001). The solid line shows the bias relation of the same functional 
form, but using the parameters a = 0.707, b = 0.35, and c = 0.80. The dotted line is the 
bias relation of Seljak & Warren (2004) calculated for a s = 0.9 and T = 0.2. The error bars, 
shown only for cr 8 = 0.6, are the error in the mean of the five realizations. Inset box: The 
galaxy bias parameter for the HOD models of Tinker et al. (2005), plotted as a function of 
(Tg. The solid lines are analytic calculations of the galaxy bias using the new bias relation. 
The dotted lines are the analytic calculations of b g using the bias relation of Seljak & Warren 
(2004). 
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Fig. 11. — The analytic HOD calculation of £ fl (r) is compared to numerical results for 
different bias and halo exclusion prescriptions. Panel (a): The calculation of £211 using the 
halo exclusion approach of Zheng (2004) is compared to N-body results. The solid line the 
calculation performed with the halo bias prescription of Sheth, Mo, & Tormen (2001). The 
dashed line is the same calculation using the modified bias in Appendix A. The error bars 
represent the error in the mean from the five realizations. Panel (b) compares the total 
calculation of £ ff (r) with the numerical results. Panels (c) and (d) compare £ 2 h and £ fl (r) 
calculated by the three new halo exclusion methods described in Appendix B to the N-body 
results. 



